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ABSTRACT 

Powerful stellar winds and supernova explosions with intense energy release in the form 
of strong shock waves can convert a sizeable part of the kinetic energy release into 
energetic particles. The starforming regions are argued as a favorable site of energetic 
particle acceleration and could be efficient sources of nonthermal emission. We present 
here a non-linear time-dependent model of particle acceleration in the vicinity of two 
closely approaching fast magnetohydrodynamic (MHD) shocks. Such MHD flows are 
expected to occur in rich young stellar cluster where a supernova is exploding in the 
vicinity of a strong stellar wind of a nearby massive star. We find that the spectrum 
of the high energy particles accelerated at the stage of two closely approaching shocks 
can be harder than that formed at a forward shock of an isolated supernova remnant. 
The presented method can be applied to model particle acceleration in a variety of 
systems with colliding MHD flows. 
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particle spectra at the maximal energies. Magnetic field 
fluctuations in the sliock vicinity may be highly amplified 
by instabilities driven by the cosmic ray c urrent and 
CR-pressure gradient in the strong sh ocks fe.g. iBell |2004| : 
iBvkov et aLll201ll : ISchure et all 120121 '). That is an impor- 
tant factor to determine the highest energies of particles 
accelerated by shocks. 

The maximum energy of accelerated particles strongly 
differs for di fferen t types of supernovae (see e.g. 
iPtuskin et all (|201(]| )). it depends on the circumstellar 
medium around the supernova progenitor star. Moreover, 
core-collapsed supernovae produced by massive stars of- 
ten occur in OB-star associations where the intense radi- 
ation of hot massive stars, powerful stellar winds and super- 
nova shocks strongly modify the interstellar environment, 
producing large hot cavities of a few tens of parsec size, 
called superbubbles. For instance, the Carina OB stars com- 
plex contains about 70 O-type stars and more than hun- 
dred B0-B3 stars c onfined in a region of about 40 pc size 
(|Gagne et al.ll201ll ). Recently, Fermi telescope detected an 
extended cocoon-shape source of gamma-ray emission as- 
sociated with a massiv e star-forming region Cygnus 0B2 
(|Ackermann et al]|201lh . The detection indicates the pres- 
ence of active particle acceleration processes in the associ- 
ation. Cygnus OB2 c ontains over 50 O-t ype stars and hun- 
dreds of B-type stars (| Wright et al.|[2O10l V Compact sources 
like binary massive stars (with the stars separated by a 
few astronomical units) may accelerate relativistic parti- 



1 INTRODUCTION 



Diffusive shock acceleration (DSA) mechanism thought 
to be operating at the fast shocks of young supernova 
remnants (SNRs) is likely responsible for the galactic 
cosmic ray (CR) accele r ation up to the "knee - region " 
energies (see e.g. iHiUad (|2005l '): lAharonian et all (|2012l ')') 
and may eve n exce ed lO^'^ eV as it was advocated by 
iPtuskin et all (|2010l ) for the case of nuclei accelerated in 
an isolated young Type lib SNRs. The basic features of 
the DSA process have b een r ev ealed in t h e pio ne ering 
paper s of lAxford et all (ll977^■ iKrvmskiil j 19771 ). lBel3 
1 19781 ') and iBlandford fc OstrikeJ (| 19781 ') where the high 
efliciency of the acceleration mechanism was shown. Effects 
of nonlinear backreaction of the accelerated particles on 
the stru cture of the superson i c shock flow were di s cusse d 
later bvlDrurv fc Volkl (1 198 ll ) IBlandford fc Eichled (Il987f) 



iBerezhko fc Krvmskiil (Il988l). iBelj (Il987l ) 



Jones fc EUisonl 
(20041), 



Blasil 



(Il99ll), iMalkov fc O'C Drurvl (l200ll) 
Amato fc Blasil (120051 ') . I'Vladimirov et all (|2008l ) and 
Reville et all (|2009l ). It has been found that the pressure 



of accelerated particles can modify the bulk plasma flow in 
the shock upstream and that may result in a substantial 
increase of the flow compression and flattening in the 



* E-mail: byk@astro.ioffe.ru 

f E-mail: peter.gladilin@gmail.com 

I E-mail: osm2004@mail.ru 



2 A.M. Bykov, P. E. Gladilm and S. M. Osipov 



cles (jEichler fc Usovl (|l993l ): iPittard fc Dougherty! (|2006l )) 
on a month time scale. Collective emission of such bina- 
ries may contribute to the gamma-ray flux observed by 
Fermi. However, if the detected gamma-ray emission is 
truly extended, it can be attributed to relativistic parti- 
cles accelera ted by multip le shocks in the superbubb le, as 



modeled by iBvkovl (1200 ll). iBvkov fc ToptvginI (|200ll ) and 
iFerrand fc MarcowithI ( 201C ) . Multiple successive interac- 



tions of particles with large shocks and rarefactions of a 
superbubble scale size occur on the time scale of about 10^ 
years. 

In the present paper we model a class of a few parsec 
size particle accelerators associated with collision of a young 
supernova shock with a fast stellar wind of a massive star. 
The modelled stage starts a few hundred years before the 
supernova shock collides with the wind termination shock 
as it is illustrated in Figure [T] At this stage the maximal en- 
ergy particles accelerated via DSA at the SNR shock reach 
the fast wind termination shock and are scattered back by 
magnetic fluctuations carried by the fast stellar wind. There- 
fore, the high energy particles that have mean free path A{p) 
larger than the distance between the two shocks L12 start 
to be accelerated by the converging fast flows. This is the 
most favorable circumstance for the efflcient Fermi acceler- 
ation. While the structure of the MHD flow in the vicinity 
of a supernova shell colliding with the stellar wind termi- 
nation shock is rather complex, it is possible to consider 
a simplified model that captures the basic features of the 
acceleration process. When the magnetic field fiuctuations 
are strong enough to provide the so-called Bohm diffusion 
regime with A(p) — ^rg{p), where p is the particle momen- 
tum, rg{p) is the particle gyroradius in the mean magnetic 
field and ^ > 1, the high energy particles bouncing between 
the converging fiows are likely to have a spectrum harder 
than that produced at an individual shock and may contain 
a sizeable fraction of the total kinetic energy of the converg- 
ing MHD fiows. Thus, we formulate a simplified nonlinear 
approach to account for the effect of the fiow modification 
by the pressure of accelerated CRs. We consider the system 
evolution stage when radii Rsh and Ravi of the two shocks 
are much larger than the distance L12 justifying a local one- 
dimensional approach. 

In § [2] we generalize for the case of two s hock fi o ws th e 
semi-analytical model orig inally developed by Blasil (|2004 ) , 



lAmato fc Blasil (|2005l ') and lCaprioli et all (|201of )" to account 
for nonlinear CR modification of a single stationary strong 
shock fiow. A two shock fiow can be stationary in the case 
of a supersonically moving star with a fast wind. In that 
case the bow shock and the termination shock are steady 
in the rest frame of the moving star. On the contrary in 
the case of a supernova shock approaching the stellar wind 
termination shock the system is non-steady. Therefore to 
model the system we constructed in § [4] a simplified time 
dependent approach valid for the shocks with fast CR accel- 
eration. If magnetic field amplification provides the efficient 
Bohm diffusion resulting in fast CR acceleration, one can 
implement CR-pressure-modified profile of the instant local 
fiow obtained in the semi-analytical solution into the time- 
dependent CR transport equation. 



2 SEMI-ANALYTIC NONLINEAR MODEL 
WITH FREE ESCAPE BOUNDARY 

Consider a model describing a population of high energy CR 
particles with A(p) > L12 in a vicinity of two approaching 
shocks with R^h ~ Rsw ^ Li2. For simplicity and the sake 
of brevity we will obtain here the solution for the shocks 
with the same parameters (velocity profile and hydrody- 
namic quantities) but this model can be easily generalized. 

In Figure [T] we illustrate a simplified scheme of the local 
fiow with two approaching shocks. At 2; < there is the 
upstream region of the shock 1 and at 2; > there is the 
upstream region of the shock 2. Shock fronts are situated at 
the a; = 0^ and a; = 0"'" positions and the CR free escape 
boundaries (FEBs) are at a; = —xq and x = xq- The CR 
FEB condition requires that 



f{-xo,p) = f[xQ,p) = 0, 



(1) 



where f{x,p) is a stationary CR distribution function. To 
derive the distribution function at the shock in the case 
of the two colliding shock fronts we employ a steady-state 
diffusion-convection equation with the CR particle injection 
rate Q(x,p) 
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(2) 



where D{x,p) is the diffusion coefficient, u{x) is the 
fluid velocity. Integrating Eq. [2] from x — —xo to x = xo 
and splitting integration region into 3 parts; from x = —xq 
to X = 0~ , from 2; = 0~ to a:: = 0"'' and from a: = 0"*" to 
a; = a::o, one can achieve the following equation 
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+ 2[4,,sc{p) - \p^^^) + Qoip) = 0, 



(3) 



where fo{p) ~ f{x = 0,p) is the distribution function at the 
shock, Qoip) — Q{0,p) is the injection rate at the shock and 
4>esc{p) = — [-D(a;,p)S^l is the escaping flux of energetic 
particles at a: = —xo and a: = a^o- 

Assuming the symmetry of the shock flow profiles and 
introducing 

1 f° 

Up = u-i— I dx{du/dx)f{x,p) , (4) 

Mp) J-^„ 

where Ui is the fiuid velocity immediately upstream (at x — 
0^ and a; = 0^), we obtain: 

The function Up is instrumental to account for the nonlinear 
modification of the fiow due to backreaction of the acceler- 
ated particles and it is i ncluded into the non-l inear calcula- 
tions (for details see e.g. lAmato fc Blasill2005l ). 

Thus the distribution function can be evaluated as 

r ^ 3(^...(p) + Qo(p)/2) 

where 
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where A(p*) = Li2). Such a simple form of the CRs distri- 
bution function can be obtained in the case of the symmet- 
ric flow. The corresponding expression for the general case 
is somewhat more complex and will be discussed elsewhere. 



3 AN APPROXIMATE SOLUTION OF 
DIFFUSION-ADVECTION EQUATION 

In this section we use an approximate solution of the 
one-dimensional diff usion-advection equation proposed by 
ICaprioli et ai] (|2010t ) with the distribution function at the 
shock fo{p) for the case of two converging shocks derived 
in the previous section. The solution was obtained by in- 
tegrating the diffusion-convection equation from —xq to an 
arbitrary point x in the shock upstream, with the FEB con- 
dition as given by Eq. ([l}. A symmetric model is considered, 
so that the solution in the upstream region 1 (a; < 0) does 
not differ from t he solution in the up stream region 2 (see 
Fig.0. Following [CaprioU et ahl (|2010l ') we used an approx- 
imation of the exact solution to the transport equation 



f{x,p) = .foexp 
uofo 



dx' 



u{x') 
D{x',p) 



W{x,p) 
Woip) 



Wo{p) ' 

where D{x,p) is the CR diffusion coefficient, 
,exp[-ij{x',p)] 



W{x,p) = ito / dx' 



D(x',p) 



V'(x,p) = 



dx' 



u{x') 



(10) 



(11) 



(12) 



(13) 



Figure 1. Simplified scheme of velocity profiles for the approach- 
ing shocks system. 



K{p,p) = exp 



"dp 3 , 
-ZT — i'^v 



1 



p u^i •■ 6 dp' 
The CR injection rate in Eq. ^ can be written as 



Q{x,p) = 



-S{p-p,„j)5(x), 



(7) 



(8) 



D{x',p) 

and Wo{p) = W{xo,p). 

These expressions are exact in the test-particle limit, 
as one can easily verify. The iterative method that has 
been used in the calculations is based on the successive ap- 
proximations to the solution f{x,p) that satisfy both the 
CR transport and the momentum-energy conservation equa- 
tions. The momentum conservation equation, normalized to 
PoUq reads 

1 



U(X) + Pcix) + P^{x) + Pg{x) = l + 



(14) 



where pi„j is particle injection momentum, 
pouo/mpUl is the gas density immediately upstream {x = 0~ 
and X = 0^), no and po are gas density and mass density at 
the shock and rj is the fraction of the particles which are con- 
sidered to be injected into the acceleration process. Hence, 
/o(p) can be presented in a simple form 



where Mo is the Mach number of the unperturbed flow. The 
normalized cosmic ray pressure 



4tt 
Spoul 



dpp'^ v{p) f{x,p), 



(15) 



fo(p) = 



Srjpo 



SirnipUpip) p^ Up{p)p^ 



(pesc P dp 



(9) 



where v{p) is the velocity of the particle. The normalized 
pressure of magnetic fluctuations ge nerated via the reson ant 
streaming instability (see Eq.(42) in lCaprioli et al.ll2009l ') 



where Up = Up/uo- 

Eq. (|9]) represents the momentum distribution function 
at the high-energy limit (A(p) > L\2) for the MHD flow 
with two identical colliding shocks with the FEBs. Note, 
that for every momentum p function fo(p) is proportional 
top~^ with correction factor Up{p). The first term in Eq. ^ 
reflects CR injection and the second term is due to the es- 
caping fiux. The second term is most important for the CR 
spectral shape at the highest energies (i.e. above p > p. 



VA l-U^jx) 

Auo C/3/2(a;) ' 



(16) 



where va = Bo/ y'A-Kpo - is the Alfven velocity, Bo - is the 
strength of the unperturbed magnetic field. Equation (|16|l 
was obtained from the stationary equation for growth and 
transport of magnetic turbulence. We used the approxima- 
tion for the turbulent energy flux F^jx) ~ 3u(x)pm(x ), as- 
suming that Va <^ u{x) (see, e.g.. ICaprioli et al ] |2009l l. The 
normalized pressure of the background gas 
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(17) 



where U{x) — u{x)/uo and 7 is the adiabatic index. 

The first iteration starts from a guess value for Ui = 
Mi/wo, which uniquely determines Pwi, Pgi and Pci via 
Eqs. ([161), GZl and (fTi)). 

To the first approximation we start with a test-particle 
guess for f{x,p), properly normalized in order to account 
for the pressure Pci, and calculate Pc{x) from Eq. 1)15^ 
and then U{x) with Eq. p4|) . Then the updated veloc- 
ity profile is used to construct a new Pm{x) and SB{x) = 
•y/87rpo^'o-PiD(a:) which is employed to update the diffusion 
coefficient D{x,p) = vpc/3eSB{x). 

According to Eq. (|10[l. a new profile of f{x,p) is ob- 
tained with the initial distribution function and the new 
U{x) and D{x,p). The procedure is iterated until a conver- 
gence is reached, i.e. until the remainders of f{x,p) and its 
normalization factor reach a preset accuracy between two 
steps. 

For an arbitrary value of (7i with the fixed list of the 
model parameters, however, the required normalization fac- 
tor can differ from 1, then the process is restarted with an- 
other choice of Ui until no further normalization is needed. 
The distribution function calculated with the value of Ui 
obtained by this approach is, by construction, the solution 
of both diffusion-convection and conservation equations for 
the two shock fiow model. 

In Fig. [2] we illustrate the CR spectrum in the limit of 
small L12 when the shocks are colliding. The proton distri- 
bution function at the shocks (dotted line) given by Eqs. ([9])- 
(|10p and the corresponding escaping flux (dashed line) given 
by Eq. (jlip were calculated with the semi-analytical ap- 
proach. The following parameters were chosen for the cal- 
culation: the shock velocity uo — 3,000 km s~^, the free 
escape boundary is at 1.0 pc on the both sides from the 
shocks, the background gas density is no = 1.0 cm~"^ 
and the Alfven velocity is va ~ 30 km s"^. Within 
the spectrum calculation we used the diffusion coefficient 
D{x,p)=5xl0'^° [p/{GeV/c)] x {5 B (x) / 100 fiG)-^ cm^ s^^ 
The ratio of the escaping and injection fiuxes is governed by 
the value of the diffusion coefficient in the upstream region 
D{x,p), and in the model it is given by the expression 



= (P) 



«o/o(p) Wo(p)' 



(18) 



At the maximal momenta of the accelerated particles the 
dimensionless ratio Eq. (|18p scales as D{xo,p). 

Note that the proton spectrum at the shock is not 
strongly modified by the cosmic ray pressure and differs 
from the case of a single modified shock simulated by 
ICaprioh et all (|201(]| ) with the same injection parameter 
£,inj ~ 3.3, corresponding to a fraction of injected particles 
»7 ~ 1.2 X 10"^ 

It is worth noting that the spectral shape of the escap- 
ing CR fiux (showed as the dashed line in Fig. [2]) differs 
from that in the case of a single strong shock. In the case 
of the approaching shocks flow the escaping particles form 
a flattened spectra at lower energies. The spectral shape is 
no longer symmetric and it departs from the parabolic law. 

In the case of a SN shock approaching a termination 
shock of the stellar wind, the system is non-steady. Never- 
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Figure 2. The CR proton spectrum fo{p) at the shock position 
(dotted line) and the escaping flux p'^4>esc{p)/uo (dashed line), 
simulated in the non-linear semi-analytical model of DSA in the 
colliding shock flow 



theless, one can make use of the non- linear steady solution 
given above to simulate the instant fiow profiles between the 
shocks, provided the CR acceleration time Ta is short com- 
pared the dynamical time scale Tdyn = L-12/usn, where its„ 
is the velocity of the SNR forward shock. The simulated in- 
stant profiles can be then implemented into the non-steady 
transport equations for CR protons and electrons to account 
for the time-dependent position of the approaching shocks. 

It is important to note that the plasma flow in between 
the converging shocks is governed by the gradient of the 
CR pressure. The CR pressure itself may be dominated by 
the highest energy part of the CR spectra. However, the 
CR protons of momenta p > propagate in between the 
shocks without scattering, and their pressure is nearly ho- 
mogeneous. In result, the flow is modified by the pressure 
gradient of the CR particles with p < p*. These have the 
acceleration time shorter than Tdyn and are mainly concen- 
trated in the vicinity of the shocks. 



4 TIME-DEPENDENT MODEL SIMULATIONS 

The non-linear modification of the colliding shock fiow by 
the accelerated cosmic ray pressure can be accounted for 
within a time-dependent model. A specific feature of the 
simulation is that the highest energy CRs of pmax S5= P ^ P* 
propagate without scattering, but their momenta are still 
nearly isotropic being scattered in the extended shock down- 
stream regions. To treat the case one can use the transport 
equation in the form of the telegraph equation. The tele- 
graph equation can be derived from the Boltzmann kinetic 
equation for a n early-isotro pic CR particle distribution (see, 
e.g., Eq.(35) in lEarl|[l974 ). The equation allows a smooth 
transition between the diffusive and the scatter-free propa- 
gation regimes. With time-dependent simulations we calcu- 
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lated the energy spectra of cosmic ray protons and electrons 
at a SNR shock approaching a strong wind of a nearby early 
type star. For relativistic electrons/positrons we account for 
the energy losses due to synchrotron and inverse Compton 
(IC) radiation. 

We solve one-dimensional transport equations for the 
pitch-angle-averaged phase space distribution function of 
protons, fp{x,p,t), and electrons, fe{x,p,t), given by 

^ (^(^'^)^) +"^P(^^^ lbeM-2y}9e] , (20) 

where pp = p'^fp, = p'^f^, y = ln{p). The dimen- 
sionless particle momentum p is expressed in the units of 
TUpC. The cooling term b(p) = —dp/dt describes th e elec- 
tron synchrotron and IC losses (see, e.g.. lKan3l201ll ). Here 
t(p) = A{p)/v is a CR particle pitch-angle scattering mean 
free time. The transport equations Eq. (1191) and Eg. (I20p are 
modi fied telegraph equations (see, e.g.. lEarll 19741 : iToptvginl 
1 19851 ). The telegraph equation embodies both the diffusive 
propagation of low energy particles p < with the mean 
free paths below L12 and the ballistic scatter free propaga- 
tion of particles with p > p*- The transport equations are 
valid for nearly isotropic particle distributions that can be 
maintained for particles with momenta pmax ^ P ^ P* that 
are scattered in the inflowing plasma (Upstream 1 and 2 
regions shown in Fig. [T]). At high energies p > p* the bal- 
listic propagation of CR particles in between the shocks is 
governed by the first terms of the transport equations. The 
diffusion coefficient in the regime saturates at Dis « cI/12. 

At the moving shocks we applied th e standard 
matching conditions used in DSA (see, e.g., iBelll 1 19781 : 
pBlandford fc Eichleil 1 19871 ') . The matching conditions are 
equating the isotropic parts of the distribution functions and 
the particle fluxes in the phase space at the shock surface 
to guarantee the particle number conservation. At the sim- 
ulation box boundaries x — ±xi, we used the free escape 
boundary conditions gii{t, —Xb,p) — and gi2{t,Xb,p) = 0. 
The simulation starts with the initial test particle distribu- 
tion function at the shocks dN{p)/dp ~ l/p'^ and dN/dt = 
at t = 0. 

The transport equations were solved with integration- 
interpolation algorithm and the standard Crank-Nicolson 
scheme. The model allows to calculate fe,p{x,p,t) at any 
position between the shocks and in the post-shock flows. 

In Fig. [3] (left panels) we present the result of calcula- 
tions of the proton distribution function. This was obtained 
in the time-dependent model simulations of the flow between 
colliding shocks. Four panels display fp{x,p,t) in the phase 
space (p, x) at four different times which correspond to the 
inter-shock distances of 0.6 pc, 0.5 pc, 0.3 pc and 0.1 pc. 
The free escape conditions are applied at the simulation box 
boundary at Xb = 0.55 pc. At each time step the flow model 
implements the modified velocity profiles, u(x,t), obtained in 
the semi-analytical iteration scheme described in paragraph 



§3. The four right panels in Fig 3 show the proton and elec- 
tron spectra, dNp{p)/dp and dNe{p)/dp, at the moving left 
shock, at the corresponding time moments, compared to the 
expected spectrum dN / {jp) / dp oc 1/p. In Fig. [3] the spectral 
evolution of CRs is clearly seen. Initially, the CR protons 
are concentrated around the shocks, and the proton spec- 
trum is close to the spectrum of the isolated SNR shock, 
dNp{p) /dp oc 1/p^. As the inter-shock distance reduces, the 
proton spectrum gets harder and the CRs concentrate in 
between the shocks. At the distance of 0.1 pc, the spectrum 
of CR protons accelerated in the two-shocks system almost 
coincides with the spectrum dN{p)/dp oc 1/p. 

The later result is consistent with the expectation that 
for the efficient acceleration process in such a system, the 
acceleration time for the highest energy particles, r^, should 
be comparable to the dynamical time Tdyn. The condition 
is satisfied if: (i) UsnL\2 > 3 • Dsn and UsnL\2 > 3 • Dsw, 
where Usn - velocity of the SNR shock. Dsn and Dsw - the 
Bohm-type diffusion coefficients for the supernova and th e 
stellar wind shocks correspondingly (see lBvkov et al.ll201lh . 
(ii) the diffusion coefficient in the interstellar medium be- 
tween the shocks Dis ^ Dsn , Dsw ■ We did not model here 
the magnetic field amplification process, but rather param- 
eterized the Bohm diffusion coefficients relying on the simu- 
lations of DS A magnetic field amplifica tion at a single shock 
performed bv IVladimirov et al.l (|20 08l). Assuming that the 
amplified magnetic field near the shock is ~ 100 fiG as it was 
inferred fro m young SNR observations recently reviewed by 
IVinkI (|2012l ). the distance between the shock fronts L12 ~ 1 
pc, we used Dis = 100 x Dsn and Dsm — Dsn for the Bohm- 
type diffusion. 



5 DISCUSSION 

The simplified non-linear model of particle acceleration in a 
converging flow between a young supernova shell and a fast 
wind of a massive early type star predicts a hard spectrum of 
CR protons confined in the fiow and CRs escaping the accel- 
erator of about a decade width in energy as shown in Fig. [S] 
In the considered case of the SNR shock and wind velocities 
of ~ 3, 000 km s~^ and about a parsec distance between the 
shocks the maximal energies accelerated protons can reach 
10^ GeV. The energy is close to the expected DSA limit 
in an isolated sup ernova shock (see e.g. iLagage fc Cesarskvl 
1 19831 : ^11^^120051 ). In individual SNRs the maximal energy 
depends strongly on the circumstellar matter of the SN pro- 
genitor and the maximal energy of the accelerated CRs can 
be reached typically at the early Sedov phase. It is impor- 
tant that if SNR is expanding in a cluster of young massive 
stars both the maximal energy and the flux of the escaping 
accelerated CRs may be non-monotonous functions of the 
SN age, and may have secondary maxima at the moments 
of close approaches of the SNR and fast stellar winds. It 
should be noted that mas sive shells that are surrou nding 
the stellar winds (see e.g. iLamers fc Cassinellil 1 19991 ) may 
be swept away by the flrst few SN explosions in a compact 
cluster of young massive stars, thus alleviating the collisions 
of supernova shells with the wind termination shocks. The 
stage of the closely approaching SNR and stellar wind flows 
that is favorable for particle acceleration typically lasts for 
300-1,000 years depending on the stellar wind velocity and 
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Figure 3. On the left: the proton distribution function f(x,p)p^ (gray scaled) as a function of the CR momentum and the position x, 
presented for four different distances between the shocks - from top to bottom: 0.6 pc, 0.5 pc, 0,3 pc, 0.1 pc. The left shock is propagating 
in the positive direction of the x- axis and the right shock is moving in the opposite direction. The shock velocities are 3,000 km . 
On the right: the proton and electron spectra at the moving left shock for the times corresponding to the four intcrshock distances given 
above. Solid line - simulated proton spectrum dNp/dp, dashed line - simulated electron spectrum dNe/dp multiplied by a factor of 100 
for the clarity of presentation. The asymptotic spectral shape dN/dp oc 1/p is shown with dotted line. 
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Figure 4. Model spectral energy distribution (SED) of the in- 
verse Compton emission from electrons accelerated at the collid- 
ing shock flow (solid line) to compare with the emission of the 
electrons accelerated by a single isolated shock of the same speed 
(dashed line). The SEDs are given in arbitrary units and calcu- 
lated for the intershock distance of 0.1 pc. The corresponding 
distribution function of the IC emitting electrons is shown in the 
bottom panel in Fig. |3] 



the termination shock radius. Therefore, such sources can 
likely contribute to the galactic CRs population. The collid- 
ing supersonic flows with two shocks also occur in the early- 
type runaway stars. Particle acceleration in a way similar to 
that discussed above may occur in between the bow-shock 
and the wind termination shock of a fast runaway early type 
star. Recently a bow-sh ock survey with 2 8 fast runaway OB 
stars was compiled by iPeri et al] (|2012h providing candi- 
date source for particle accelerators. The velocities of run- 
away early type stars are well below that of the supernova 
shocks discussed above and therefore the maximal energies 
of accelerated particles are expected to be lower than that 
in SNR-SW collision. However the fast runaway early type 
stars may still be sources of non-thermal emission with hard 
photon spectra. 

The simplifled non-linear model discussed in the pa- 
per predicts a hard high energy end of CR particle spec- 
trum. The simulated spectrum of electrons is shown in Fig. [3] 
(dashed curve). The shape of the electron spectrum is sim- 
ilar to that of protons, but the synchrotron/IC cooling ef- 
fects result in the spectrum cut-off at lower energies. The 
electromagnetic radiation produced by hard spectra of elec- 
tron/positron confined in the accelerator can be potentially 
observed, thus constraining the physical parameters of the 
observed objects. In Fig. [4] we present a model spectrum 
of the inverse Compton (IC) emission from the electrons 



accelerated in a colliding shock flow comparing to the IC 
spectrum of a single shock flow mod e led earlier (see e.i 
iGaisser et al] 119981 : [Baring et al.lll999l : iBvkov et all I200C 
To construct the IC spectrum we used the same emission 
model which was applied to model high energy emission 
of the supernova remnant IC 443 described in detail in 
iBvkov et al] (|2000l ). SNR IC 443 is likely located in Gem 
OBI - young massive star association. Therefore in addition 
to the local interstel lar photon spectrum which following 
iMathis et all (| 1983h was approximated by the sum of di- 
luted blackbody spectra with temperatures 2.7 K, 7500 K, 
4000 K and 3000 K a major infrared component resulting 
from the dust heated in OB association must be accounted 
for. In t he vicinity o f IC 4 43 according to IRAS observa- 
tions of ISaken et al.l (| 19921 ') we approximated the infrared 
component with two temperatures flt (185 K and 34.3 K 
components). If the shock collision region is located at the 
distance further than about 2 parsecs from the young mas- 
sive star then the UV emission of the star does not dominate 
the IC losses of relativistic electrons. 



The electron distribution used in the IC spectrum com- 
putation is shown in the bottom panel in Fig. [3] The 
electron/positron distribution function fe{x,p,t) was com- 
puted with the time-dependent model Ea. H20p that accounts 
for non-linear flow modification in the shocks vicinity as 
discussed in the previous section, and the electron syn- 
chrotron/IC losses included. It is clearly seen that due to 
the hard spectrum of accelerated electrons up to 10* GeV 
the emissivity of the SNR-SW system is higher than that 
for the single SNR shock. TeV-sources with very hard spec- 
tra are not always easily identified with their GeV coun- 
terparts and some of these may comprise the population of 
so-called "dark accelerators". The prototype of the dark ac- 
celerators was TeV J20 32-I-4130 discovered a decade ago by 
lAh aronian et al.1 (|2002l ). though the source is possibly asso- 
ciated with the pulsar 2FGL J2032.2-f 4126 detected later 
in GeV regime by Fermi observa tory. Another high energ; 
source was found by H.E.S.S. (see lAbramowski et al.ll2012V ' 
in the vicinity of the young massive stellar cluster Wester- 
lund 1. The sources originating in the colliding supersonic 
flows considered above may contribute to the VH E emis- 
sion of young massive s tar associations (see, e.g., iBvkovl 
I2OOII : [Torres et al] l2004h and the star-burst galaxies that 
are also known to be bri ght gamma-ray sources (see, e.g., 
lAbramowski et al.|[2012bl ). The main goal of the paper is to 
discuss the unique features of particle accelerators associated 
with the colliding shock flows. To model the broadband non- 
thermal emission spectra of the sources one has to discuss in 
detail matter and magnetic fleld distributions in the complex 
flows. The structure of magnetic fleld in the w inds of young 
massiv e stars is under study, see for a review IWalder et al.l 
(I2OI2I ). In the IC spectra modeling we assumed that the CR 
electrons are the test particles that are propagating in the 
flow modified by the CR protons. We did not consider here 
any specific electron injection model to calculate their abso- 
lute fiuxes and therefore the IC spectra in Fig. U are shown 
in arbitrary units. We will present a more detailed emission 
spectra modeling of colliding shock flows elsewhere. 
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APPENDIX A: MOMENTUM-ENERGY FLUX 
CONSERVATION 

In this section we check the accuracy of the quasi-steady 
non-linear model used, since it is relying on an the ap- 
proximation of the CR distribution function as given by 
Eq. (|10|) . In a steady-state system, mass, momentum and 
energy fluxes must be constant in space obeying the follow- 
ing conservation laws 



p(x)u(x) = pouo, 
$p(x-) = poul + P< 



sO, 



$e(x) 



(Al) 
(A2) 
(A3) 



where p and u are the mass density and the flow velocity, 
$p(a;) and $e{x) are the fluxes of the s-component of the 
momentum and energy in the x-direction. The index "0" 
marks the far upstream values. 

The CR energy and momentum fluxes are calculated 
from the CR distribution function. In the paper we assumed 
that the CR energy and momentum are dominated by pro- 
tons, and the CR electrons are treated as test particles that 
are propagating in the flow modified by CR protons. A 
weakly anisotropic CR distribution function of the accel- 
erated particles can be approximated as 



/(r,p) 



1 

4% 



N{r,p) + ^vJ(r,p) 



(A4) 



where N{r,p) is the isotropic part of the distribution func- 
tion, J(r,p) is the particles flux (e.g. lToptve:in|[T985l ). The 
2;-component of the flux can be written as: 



J:c{x,p) = -D{x,p) 



ON pdN 



u(x) . 



dx 3 dp 
The energy flux "l>cr(a;) is deflned as 



(A5) 
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where Pm{x) — poUgP-ujix) is the pressure of the magnetic 
turbulence g enerated via resonant streaming instability (see 
Eq. (HSl) and lCaprioIi et alj (|2009h l. 

Therefore, the energy conservation law for the steady 
non-linear model can be written as 



1 3 '~y 

-p(x)u {x)+3u{x)pu,{x)-\ -u(x)pg{x) 

2 7 — i 



-^crix) 



(All) 

const{A12) 



where ^cr{x) is calculated from Eq. (|A9|) . 

Pg{x) = poulPg{x) 

is the pressure of the background gas and 



E ■ 4>esc{p)dp 



is the energy flux of the escaping CRs. 

In Figure lAll we illustrate the accuracy of the energy 
conservation in the employed model of particle accelera- 
tion. The dimensionless spatial coordinate is normalized to 
the FEB position distance xo- The energy flux conserva- 
tion accuracy is about 1%. Some small inconstancy of the 
full energy flux can be attributed to the particular choice of 
the approximation of the CR distribution function given by 
Eq. (|1U|) that was used in the non-linear model. 



Figure Al. The energy conservation check. The shock velocity 
profile U (x) is given for the fiow of Mach number Mq = 100. 
The shock is located at a; = and the position of the free es- 
cape boundary is at a; = 1.0. The cosmic ray pressure Pc{x), the 
pressure of the background gas Pg{x) and the pressure of mag- 
netic fluctuations Pw{x) are normalized to PqUq. The energy flux 
^e{x) is normalized to pou^. 



^,r{x)^ j K{p)v4p)f{x,p)d'p, (A6) 

where K(p) = E — m(? is the kinetic energy of a particle, 
E = (pc)^ -I- m^c* is the full energy of the particle, v^, = 
V ■ cosO, 9 is the angle between the particle momentum and 
the X-axis, and v = 

Substituting the a;-component of the distribution func- 



tion Eq. (|A4|) into the equation Eq. (|A6|I . one can obtain 



^cr{x) 



N{x,p) -f -cos6Jx{x, 



K{p)v cose ^.{A7) 



Integrating over the angle one can get the CR energy flux 



<E.c.(x-)= / K{p) 



ON pdN ' 



p'^dp. (A8) 



Then, flnally, the expression for the CR kinetic energy flux 
is 



$cr{x) 



p^dp 



-Dix,p)^K{p) + {^+ K{p))Nu 



.(A9) 



The energy flux of the magnetic turbulence is expressed by 
^^{x):^3u{x)p^{x), (AlO) 



